% n=2, two I(1)


function [A, r, p, x]=simul2_3(T)

% first I(1) component

w=randn;

for t=2:T,
    w(t)=w(t-1)+randn;
end


% second I(1) component

u=randn(1,2);

for t=2:T,
    u(t)=u(t-1)+randn;
end



n=2;

Q=ort(2);
x=Q*[u;w];

pmax=6;
BIC=zeros(1,pmax);
    
for p=1:pmax,
    
    [phi_{p}, C_{p}]=reg(x,p);
    
    BIC(p)=log(det(C_{p}))+p*n^2*log(T)/T;
    
end

[mB, p]=min(BIC);
C=C_{p};
phi=phi_{p};


for i=1:n
    for j=1:n
        
        if i==j,
            coef=[1 -phi(i+(0:p-1)*n,j)'];
        else
            coef=[0 -phi(i+(0:p-1)*n,j)'];
        end
        A{j,i}=coef(end:-1:1);


    end
end

F=[phi';eye((p-1)*n) zeros((p-1)*n,n)];

r=eig(F);
    
    
    
